Пространственно временные модели соотяний. Фильтр Калмана

Эффективная и очень популярная в приложениях модель, оценивающая вектор состояний динамической системы, используя ряд неполных и зашумленных измерений была предлдожена Рудольфом Калманом в 1960 году. Модель получила название фильтр Кальмана.Фильтр Калмана — это, наверное, самый популярный алгоритм фильтрации, используемый во многих областях науки и техники. Благодаря своей простоте и эффективности его можно встретить в GPS-приемниках, обработчиках показаний датчиков, при реализации систем управления в системах автопилота, управления движения подводных лодок, в финансовой сфере и т.д. Пусть \(y_t\) - n-мерный вектор наблюдаемый в момент t. Пусть динамика модели описывается r-мерным вектором \(\xi_t\) (возможно не наблюдаемый), который носит название вектора состояний (state vector). Пусть динамика процесса \(y_t\) описывается следующей систмой уравнений:\[\xi_{t+1}=F\xi_t+v_t\] \[y_t=A'x_t+H'\xi_t+w_t\] где \(F,A',H'\) матрицы параметров модели размера #(rr),(nk), (nr)$ соответственно, \(x_t\)- k-мерный вектор экзогенненных (предопределенных) переменных. Первое из уравнений носит название “уравнение состояний”(State equation),второе “уравнение наблюдений”, а про вектора \(w_t,v_t\)известно, что их ковариационные матрицы \(Q,R\) соответственно и для любых моментов времени \(t,\tau\) \[cov(v_t,w'_{\tau})= 0\] Системы уравнений обычно используются для описания конечного числа наблюдений \(y_1,...,t_T\) для которых начальное значение вектора состояний \(\xi_1\) задано. Также будем предполагать, что \[E(v_t,\xi_1)=0; t = 1,2,...,T\],\[E(w_t,\xi_1)=0; t = 1,2,...,T (*)\]

Из первого уравнения системы следует, что \[\xi_t=v_t+Fv_t{-1}+F^2v_{t-2}+...+F^{t-2}v_{2}+F^{t-1}\xi_1\] \[t=2,3,...,T\] Из прдположений (*) следует, что \[E(v_t,\xi'_{\tau})=0; \tau = t-1,t-2,...,1\] Аналогично из предположений модели немедленно получаем: \[E(w_t,\xi'_{\tau})= 0; \tau = 1,2,...,T\] \[E(w_t,y'_{\tau})= E(w_t,(A'x_{\tau}+H'\xi_{\tau}+w_t)')=0; \tau = t-1,t-2,...,1\] \[E(v_t,y'_{\tau})= 0; \tau = t-1,t-2,...,1\] Для примера рассмотрим одномерный AR(p) процесс.
Пусть \[y_{t+1}-\mu= \phi_1(y_{t}-\mu)+...\phi_p (y_{t-p+1}-\mu)+\epsilon_t\] Тогда это можно переписать в рамках рассматриваемой модели.

Уравнение состяний будет:

\(\begin{equation*} \left[ \begin{matrix} y_{t+1}-\mu \\ y_{t}-\mu\\ ....\\....\\ y_{t-p+2}-\mu\end{matrix} \right] =\left[ \begin{matrix}\phi_1&\phi_2&....&\phi_p\\ 1&0& ...&0\\0&1&...&0\\.....\\0&...&1&0 \end{matrix} \right]\left[ \begin{matrix}y_{t}-\mu\\y_{t-1}-\mu\\....\\y_{t-p+1}-\mu \end{matrix} \right] \end{equation*}\)

Уравнение наблюдений

\[\begin{equation*} y_t=\mu+\left[ \begin{matrix} 1&0&0&...&0\end{matrix} \right] \left[ \begin{matrix} y_t-\mu\\y_{t-1}-\mu\\...\\y_{t-p+1}-\mu \end{matrix} \right]\end{equation*}\] Итак \[\begin{equation*} \xi_t=\left[ \begin{matrix}y_t-\mu\\y_{t-1}-\mu\\...\\ y_{t-p+1}-\mu\end{matrix} \right]\end{equation*}\] \[\begin{equation*}F=\left[ \begin{matrix}\phi_1&\phi_2&...&\phi_p \\1&0&...&0\\0&1&...&0\\0&...&1&0\end{matrix} \right]\end{equation*}\] \[\begin{equation*}v_{t+1}=\left[ \begin{matrix}\epsilon_{t+1}\\0\\..\\0\end{matrix} \right]\end{equation*}\] \[\begin{equation*}Q=\left[ \begin{matrix}\sigma^2&0&...&0\\0&0&..&0\\ ...\\0&0&...&0\end{matrix} \right]\end{equation*}\] \[A'=\mu,x_t=1\] \[\begin{equation*}H'=\left[ \begin{matrix} 1&0&...&0 \end{matrix} \right]\end{equation*}, w_t= 0,R=0\]

Все параметры модели заданы

В качестве второго примера рассмотрим модель скользящего среднего первого порядка. \[y_t = \mu+\epsilon_{t}+\theta \epsilon_{t-1}\]

Уравнение состояний

\[\begin{equation*}\left[ \begin{matrix} \epsilon_{t+1}\\ \epsilon_t\end{matrix} \right]=\left[ \begin{matrix}0&0\\1&0 \end{matrix}\right]\left[ \begin{matrix} \epsilon_{t}\\ \epsilon_{t-1}\end{matrix} \right]+ \left[ \begin{matrix} \epsilon_{t+1}\\ 0\end{matrix} \right]\end{equation*}\]

Уравнение наблюдений

\[y_t=\mu+\begin{equation*}\left[ \begin{matrix} 1&\theta\end{matrix} \right] \left[ \begin{matrix} \epsilon_t\\\epsilon_{t-1}\end{matrix} \right]\end{equation*}\]

Таким образом

\[\begin{equation*} \xi_t=\left[ \begin{matrix}\epsilon_t\\ \epsilon_{t-1} \end{matrix} \right]\end{equation*}\] \[\begin{equation*}F=\left[ \begin{matrix}0&0\\1&0\end{matrix} \right]\end{equation*}\]

\[\begin{equation*}v_{t+1}=\left[ \begin{matrix}\epsilon_{t+1}\\0\end{matrix} \right]\end{equation*}\]

\[\begin{equation*}Q=\left[ \begin{matrix}\sigma^2&0\\0&0\end{matrix} \right]\end{equation*}\] \[A'=\mu,x_t=1\] \[\begin{equation*}H'=\left[ \begin{matrix} 1&\theta\end{matrix} \right]\end{equation*}, w_t= 0,R=0\]

Модель ARMA(p,q)

Представление модели произвольного порядка \(p,q\). Пусть \(r = max[p,q+1]\) Недостающие до \(r\) параметры \(\phi_k,\theta_k\) доопределим нулями.

\[y_{t+1}-\mu= \phi_1(y_{t}-\mu)+...\phi_r (y_{t-r+1}-\mu)+\epsilon_t+\theta_1\epsilon_{t-1}+...\theta_{r-1}\epsilon_{t-r+1}\]

Уравнение состяний будет:

\[\begin{equation*} \xi_{t+1} =\left[\begin{matrix}\phi_1&\phi_2&....&\phi_p\\ 1&0& ...&0\\0&1&...&0\\.....\\0&...&1&0 \end{matrix} \right]\xi_t +\left[ \begin{matrix}\epsilon_{t+1}\\0\\....\\0 \end{matrix} \right] \end{equation*}\]

Уравнение наблюдений

\[y_t=\mu+\begin{equation*}\left[ \begin{matrix} 1&\theta_1&\theta_2&...&\theta_r&\end{matrix} \right]\xi_t\end{equation*}\]

Вывод фильтра Калмана

Повторим определение модели. \[\xi_{t+1}= F\xi_{t}+v_{t+1}\] \[y_{t}= A'x_t+H'\xi_t+w_t\]

\(cov(w_t,w_{\tau})= R\) при \(t=\tau\) иначе 0.

\(cov(v_t,v_{\tau})= Q\) при \(t=\tau\) иначе 0.

Пусть нам уже доступны данные \(x_1,x_2,....,x_T\) и \(y_1,y_2,...,y_T\). Также допустим нам уже известны оценки для \(F,A',H',Q,R\) позднее рассмотрим способ получения этих оценок.

Пусть \[\tilde{\xi}_{t+1|t}=E[\xi_{t+1}|Y_t]\]
где \(Y_t= x_1,x_2,...x_t;y_1,y_2,...,y_t\) линейный прогноз состояния на момент момент времени \(t+1\) в момент \(t\), при условии, что \(Y_t\) уже известно. Фильтр Калмана получает эти прогнозы последовательно \(\tilde{\xi}_{1|0},\tilde{\xi}_{2|1},...,\tilde{\xi}_{T|T-1}\) сопоставляя каждому из этих прогнозов матрицу минимума квадрата ошибок(MSE). \[P_{t+1|t}=E[\xi_{t+1}-\tilde{\xi}_{t+1|t})(\xi_{t+1}-\tilde{\xi}_{t+1|t})']\]

Рекурсия начинается с \(\tilde{\xi}_{1|0}\), которая в начале есть просто \[\tilde{\xi}_{1|0}=E(\xi_1)\] сопоставим ей MSE матрицу \[P_{1|0}=E[(\xi_1-E\xi_1)(\xi_1-E\xi_1)']\]

Для примера в MA(1) у нас было \[\begin{equation*} \xi_1=E\left[ \begin{matrix}\epsilon_1\\ \epsilon_{0}\end{matrix} \right]=\left[ \begin{matrix}0\\ 0\end{matrix} \right]\end{equation*}\] поэтому \[\begin{equation*}P_{1|0}=E\left[ \begin{matrix}\epsilon_1\\ \epsilon_0\end{matrix}\right] \left[ \begin{matrix}\epsilon_1&\epsilon_0\end{matrix}\right]=\left[ \begin{matrix}\sigma^2&0\\0&\sigma^2\end{matrix}\right]\end{equation*}\]

В общем случае, при предположении, что процесс \(\xi_t\)-cтационарный \[E(\xi_{t+1})= FE(\xi_t)\]. Переносим в левую часть, и последнее соотношение перепишем \[(I-F)E(\xi_t)=0\] Здесь мы воспользовались тем, что для стационарного процесса \(E(\xi_{t+1})= E(\xi_t)\). Уравнение имеет единственное решение если \(I-F\) не вырожденная матрица. Решение это \(E(\xi_t)= 0\)

Безусловная дисперсия вектора \(\xi\) \[E(\xi_{t+1}\xi'{t+1})= E(F\xi_t+v_{t+1})(F'\xi'_t+v'_{t+1})=FE(\xi_t\xi'_t)F'+E(v_{t+1}v'_{t+1})\] Положим \(\Sigma=E(\xi_t\xi'_t)\),тогда \[\Sigma=F\Sigma F'+Q\]

Решение которого \[vec(P_{1|0})=vec(\Sigma)= [I-FF']^{-1}vec(Q)\] Итак вычисление \(P_{1|0}\) на первой итерации закончено.

Прогнозирование \(y_t\)

Мы получили \(\tilde\xi_{1|0}\) и \(P_{1|0}\) далее последовательно получим по индукции \(\tilde\xi_{t|t-1}\) и \(P_{t|t-1}\) цель вычислить \(\tilde\xi_{t+1|t}\) и \(P_{t+1|t}\)

Вычисляем прогноз \[\tilde y_{t|t-1}=E(y_t|x_t,Y_{t-1})\] \[\tilde y_{t|t-1}=E(y_t|x_t,Y_{t-1})=A'x_t+H'\xi_t=A'x_t+H'\tilde\xi_{t|t-1}\] Ошибка прогноза будет \[E(y_t-\tilde y_{t|t-1})(y_t-\tilde y_{t|t-1})'=E[H'(\xi_t-\tilde \xi_{t|t-1})(\xi_t-\tilde \xi_{t|t-1})'H]+E[w_t w_t']\]

Здесь мы воспользовались тем, что члены при перекрестном умножении вида \[E(\xi_t-\tilde \xi_{t|t-1})w_t=0\]

В итоге получим, что дисперсия ошибки прогноза \[E(y_t-\tilde y_{t|t-1})(y_t-\tilde y_{t|t-1})'=H'P_{t|t-1}H+R\]

Data <- read.csv("c:/tmp/brent.csv",sep=",")
dt<-  Data[,2]
matplot(dt,type="l",col = "blue",main = "rub/usd",lwd=2)

Ряд не стационарен, перейдем к разностям

dtdif <- diff(dt)

matplot(dtdif,type="l",col = "blue",main = "Diffrenced rub/usd",lwd=2)

Посмотрим ACF и PACF

acf(dtdif,col="blue",lwd  =5)

pacf(dtdif,col="blue",lwd  =5)

PACF идентифицирует модель как AR(2) подгоним ее и представим в виде уравнений состояния и наблюдений фильтра Калмана Уравнение состояний \[\xi_{t+1} = F\xi_t+u_t\] уравнение наблюдений \[y_t=H\xi_t+v_t\]

fit2 <- arima(dtdif, c(2, 0, 0))

Матрица переходоа\(F\)

fit2$model$T
##            [,1] [,2]
## [1,] -0.1566701    1
## [2,] -0.1245452    0

Матрица \(H\)

fit2$model$Z
## [1] 1 0

Текущее состояние \(\tilde\xi_{T-1:T}\)

fit2$model$a
## [1]  1.794680e-01 -5.828181e-05

Сглаживание при помощи фильтра Калмана

klike<- KalmanSmooth(dtdif,mod=fit2$model)
tot<-cbind(klike$smooth[,1],dtdif)
matplot(tot,type ="l",lwd=2,main = "Smoothing")

Прогнозирование \(y_t\) на \(l=10\) шагов вперед

forec <- KalmanForecast(n.ahead = 10L, mod=fit2$model, update = FALSE)
forec
## $pred
##  [1] -2.817554e-02 -1.793761e-02  6.319414e-03  1.243980e-03 -9.819472e-04
##  [6] -1.090053e-06  1.224676e-04 -1.905124e-05 -1.226799e-05  4.294768e-06
## 
## $var
##  [1] 1.000000 1.024546 1.034545 1.035783 1.035831 1.035861 1.035861
##  [8] 1.035862 1.035862 1.035862

Оценка параметров

При предположении о распределении \(v_t,w_t\) как многомерном нормальном с нулевым математическим ожиданием и ковариационными матрицами \(Q,R\) соответственно условное распределение \(y_t\) при условии, что известно \(x_,Y_{t-1}= (x_1,x_2,...x_{t-1};y_1,y_2,...,y_{t-1}\) будет нормальным

\[y_t|x_t,Y_{t-1}\sim N(A'x_t+H'\xi_{t|t-1},(H'P_{t|t-1}H+R))\] поэтому \[f(y_t|x_t,Y_{t-1})=(2\pi)^{n/2}|H'P_{t|t-1}H|^{-1/2}*exp(-1/2(y_t-A'x_t-H'\xi_{t|t-1})'(H'P_{t|t-1}H+R)^{-1}(y_t-A'x_t-H'\xi_{t|t-1}); t = 1,2,...T\] Отсюда приходим к логарифму функции правдоподобия \[\sum_{t=1}^Tln(f(y_t|x_t,Y_{t-1})\], которое максимизируется численно относительно параметров модели \(F,Q,A,H,R\)

Фильтр Калмана пытается оценить состояние системы в текущий момент времени, основываясь на своей оценке в предыдущий момент времени, а также на новых результатах измерений. Обозначим априорную оценку (до получения результатов измерения y k